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SUMMARY 

This  paper  presents  a  method  Jor  decomposing  bulk  stress  data  into  individual  stress  components. 
This  has  applications  in  the  field  of  thermoelastic  stress  analysis  where  the  raw  data  are  related  only  to 
the  sum  of  the  principal  stresses  f referred  to  as  bulk  stress).  By  considering  equilibrium,  and  with  some 
knowledye  of  the  form  of  the  solution  and  boundary  conditions,  it  is  shown  that  bulk  stress  data  for  a  two 
dimensional  body  may  be  separated  into  stress  components  by  mean  of  a  least  squares  method. 
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1.  INTRODUCTION 

Whilst  the  thennoelastic  effect  has  been  known  for  well  over  a  century,  it  is  only 
within  recent  years  that  this  phenomenon  has  been  exploited  as  a  means  of  ex¬ 
perimental  stress  analysis.  Equipped  with  an  infrared  sensor,  together  with  some 
precision  scanning  mechanisms  and  clever  signal  processing  hardware,  a  system 
known  as  SPATE  (Stress  Pattern  Analysis  by  measurement  of  Thermal  Emis¬ 
sion),  can  form  a  raster  picture  of  the  temperature  fluctuation  of  body  under 
cyclic  loads.  Since  such  thermal  data  are  directly  related  to  the  changes  in  bulk 
stress  ( cXj-jt  -I-  +  <r22),  this  system  offers  a  powerful  means  of  obtaining  full- 

field  non-contact  measurements  of  stresses  within  structural  components.  As  a 
result,  the  use  of  SPATE  has  become  increasingly  widespread,  and  many  success¬ 
ful  applications  may  be  found  in  the  literature  (e.g.,  [l]-[3] ) .  However,  because 
temperature  and  bulk  stress  are  scalar  values,  the  vectorial  nature  of  stress  is  lost 
from  data  obtained  by  such  techniques.  Whilst  the  bulk  stress  information  may 
be  useful  for  many  qualitative  inspection  purposes,  quantitative  evaluations  are 
possible  only  for  cases  where  it  is  known  that  the  stress  field  is  dominant  in  one 
direction.  Taking  an  extreme  case  as  an  example,  it  may  be  shown  that  SPATE 
would  not  be  able  to  distinguish  an  unloaded  body  and  one  which  is  subjected  to 
pure  shear.  The  failure  of  SPATE  to  respond  to  pure  shear  was  clearly  shown  in 
the  experiment  of  Stanley  and  Chan  [*!]. 

Because  the  knowledge  of  individual  stress  components  is  important  in  many 
practical  situations,  much  work  has  been  done  in  the  fields  of  holographic  interfer¬ 
ometry  to  separate  principal  stress  components  from  isopachic  and  isochromatic 
fringe  data  (e.g.  [5-7] ) .  In  thermoelastic  stress  analysis  however,  only  the  isopachic 
data  may  be  obtained  which,  by  themselves,  do  not  provide  enough  information 
for  determining  the  stress  components.  On  the  other  hand,  it  is  not  true  that 
the  selection  of  stress  values  which  satisfy  a  given  set  of  isopachic  data  can  be 
totally  arbitrary.  Phis  is  because  the  permissible  stresses  must  satisfy  equilibrium 
as  well  as  known  boundary  conditions.  It  is  shown  in  this  paper  that,  at  least 
for  two-dimensional  problems,  the  imposition  or  these  conditions  can  be  used  to 
determine  the  stress  components  associated  with  a  given  set  of  isopachic  data. 


2.  THEORY 

For  a  two-dimensional  isotropic  material  under  elastic  deformation  in  the  absence 
of  any  body  force,  it  is  well  known  that  the  solution  for  stresses  may  be  expressed 
in  terms  of  a  potential  <t>,  such  that 

d2*  d 2#  d 2$>  ,nt. 

a*x~  dy >’  ***  a*2’  <Tr"  “  (21) 

where  <rxx,  <rn  and  <Tzy  are  the  normal  and  shear  stress  components,  and  4>  must 
satisfy  the  bi-harmonic  equation 

a4*  a4*  a4*  „  , 

ax4  +  2ax3ay2  +  a^  “  ‘  (22) 

For  a  simply  connected  body,  the  specification  of  two  conditions  of  4>  around 
the  boundary  forms  a  well  posed  problem,  and  it  may  be  shown  that  the  necessary 
boundary  conditions  for  4>  may  be  determined  from  two  known  components  of 
stress  at  each  point  along  the  boundary.  In  such  a  case,  the  stress  field  may  be 
uniquely  determined,  and  there  is  thus  no  need  for  additional  data.  However, 
accurate  boundary  conditions  may  be  unknown  in  many  real  situations,  and  the 
analyst  is  forced  to  resort  to  simplifications  which  may  or  may  not  be  realistic. 
The  assumption  of  a  uniformly  distributed  stress  profile  across  the  boundary  of  a 
component  under  consideration  is  one  such  example.  With  the  availability  of  other 
stress  related  information,  such  as  experimentally  obtained  bulk  stress  data,  the 
requirement  for  a  complete  set  of  the  boundary  conditions  may  be  relaxed.  The 
problem  now  becomes  one  which  involves  finding  a  solution  to  eqn  (2.2)  subject 
to  certain  known  boundary  conditions  (e  g.,  making  use  of  free  boundaries)  and 
its  nearness  to  I  lie  bulk  stress  data. 

As  an  illustration,  consider  a  rectangular  segment  of  a  specimen  hounded  by 
0  <  x  <  (,  and  —0.5  <  y  <  0.5.  bet  the  sides  y  =  ±0.5  be  free  edges  and  the 
remaining  two  boundaries  be  loaded  by  normal  stresses  which  may  be  functions  of 
y.  Without  loss  of  generality,  considering  only  the  symmetrical  (about  the  x-axis) 
case,  the  analytical  solution  for  <J>  may  be  found  in  Timoshenko  and  tJoodier  [8]. 


<t>  =  at-  +ur_2,-r 


-•)  tan  -  cos '>■)!/  +  2-jysin  2 ~,y) . 


in  which  o0  ami  a  are  real  coefficients,  and  *,  satisfies  the  equation 

sin  2*)  +  2y  =  0. 


The  non-trivial  roots  of  e<jn  (2.4)  are  complex,  and  appear  in  conjugate  pairs. 
Also,  if  7  is  a  root,  so  is  —7.  For  simplicity,  let  us  further  assume  that  the  stresses 
are  known  to  be  uniformly  distributed  as  1  — *  00  so  that  only  values  of  7  which 
have  positive  real  parts  need  be  considered.  For  the  stresses  to  be  real,  the  general 
expression  for  $  may  be  written  as 


*  =  a»4"  +  (— 7»  tain*  cos27tj  +  27tj/sin27*y) 

*=l 

rx. 

+  '%2akc--'’>T  (- 7t  tan  7* cos 2^ky  +  2^kys\n  2 7*St). 

k=  1 


(2.5) 


where  ak  and  7t  are  the  conjugates  of  ak  and  -)k  respectively. 

The  expressions  for  the  stress  components  may  be  obtained  by  direct  differ¬ 
entiation  of  e<pi  (2.5)  and  it  may  be  shown  that  the  bulk  stress  S  is  given  by 


S  =  a  Tr  +  (r„y 

■x  X  0) 

=  n0  cos 2-) A  //  +  <*<>s2 

*=i  fr=i 

In  general,  the  coefficients  <ik  would  be  dependent  on  the  distribution  of  ir,..,. 
and  cr,,,  at  j  =  0.  However,  suppose  that  this  stress  distribution  is  not  known 
a  priori,  but  experimental  data  on  bulk  stresses,  such  as  those  obtained  from 
SPATE,  are  available.  The  problem  of  determining  individual  stress  components 
for  this  example  then  essentially  involves  the  determination  of  the  coefficients  ak 
(and  also  the  number  of  fit's  to  sufficiently  represent  tile  solution)  in  <><|n  (2.5) 
subject  to  some  best-lit  criterion  such  that  the  difference  between  .S  and  the  array 
of  experimentally  observed  bulk  stresses  S1  is  minimised. 
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3.  TEST  CASES 

To  test  the  viability  of  the  proposed  scheme,  the  stress  field  was  first  solved  directly 
for  the  case  where  a  known  distribution  of  stresses  at  the  boundary  x  =  0  is  applied. 
Th<*  bulk  stresses  at  discrete  points  on  a  A'x  x  Nt  grid  were  then  generated  and, 
with  white  noise  of  various  degree  added,  served  as  a  simulated  SPATE  output. 
This  set  of  data  was  then  used  in  the  inverse  problem  where  the  coefficients  a* 
were  estimated  and  the  resulting  stress  components  compared  to  those  obtained 
in  the  direct  problem. 

3.1  Solving  for  y 

In  order  to  solve  either  the  direct  or  inverse  problem,  the  solution  to  otjn  ( 2.4 ) 
must  firstly  be  sought.  Because  of  its  nun-algebraic  nature,  a  numerical  algorithm 
using  a  Newton- Raphson  scheme  was  used.  The  asymptotic  form 

..  /3?r  \  ,  lu(d<-  -f  ;l)?r 

(lmi  =  I  —  4- for  J  4- 1— - - -  (3.1) 

was  lists  I  to  generate  inil  ial  estimates  for  t  he  iterative  solut  ion  procedure.  This  and 
all  subsequent  computations  were  done  on  an  EI.XSI  (MOO  and  double  precision 
programming  was  used  throughout . 


3.2  The  Direct  Problem 


The  first  step  in  the  direct  problem  was  to  decide  on  a  stress  distribution  at  the 
boundary  x  =  0.  Serving  as  a  relatively  severe  test  case,  a  step  function  of  the 
form 


<7^=100  -  1/3  <  if  <1/3 

<t  j-j.  =  0  elsewhere 


(3.2) 


was  chosen,  where  the  normal  stresses  <rzr  at  t  =  0  may  be  expressed  as 


N 

<rxx  =  tan  -)*  +  )cos2->*y  -  8->2ysin  2>*y} 

jfc=] 


N 

+  ian  ^  +  8^)cos  2Tj “  *T?ysin  2T ky}y 

k  =  l 


(3.3) 


where  N  is  sufficiently  large  to  represent  the  function  expressed  by  eqns  (3.2). 


The  coefficients  a*  were  then  solved  by  matching  the  first  N  Fourier  com¬ 
ponents  of  both  sides  of  eqn  (3.3).  It  was  found  that  N  =  15  was  adequate  to 
represent  the  required  step  function. 

A  regular  grid  system  of  241  x  121  points  was  used  to  cover  the  region  0  < 
x  <  2  and  the  bulk  stress  at  each  grid  point  was  generated  using  the  truncated 
form  of  eqn  (2.6).  Random  errors  (up  to  ±5  units)  were  then  added  to  the  data  to 
form  a  simulated  SPATE  output  SE .  Contours  of  the  generated  bulk  stress  and 
the  simulated  SPATE  output  SE  (over  the  region  0  <  x  <  1)  are  shown  in  Fig.  1. 


3.3  The  Inverse  Problem 

As  ment  ioned  earlier,  the  inverse  problem  involves  the  esl  imatiou  of  I  he  coefficients 
nk  given  a  set  of  bulk  stress  data  SE .  Phis  at  first  appears  straightforward  and 
a  least  squares  approach  seemed  appropriate  for  handling  the  noise.  However, 
problems  can  arise  due  to  under-parameterisation  (the  under-estimation  of  the 
number  of  terms  necessary  to  represent  the  solution)  as  well  as  poor  conditioning 
of  the  least  squares  matrix.  These  problems,  ami  the  suggested  remedies,  are 
discussed  in  the  following. 

The  bulk  stress  may  be  rewritten  in  the  form 


.S'(x,  y )  =  a,  +  621-1  ros2")‘  l/  +  cos 2 ^ky) 

*=  1 

.v 

+  ‘Ylb7k  {e~2~’lT  Ity  ~  c~2'*T  ™s2*iky) . 

*=t 


(3.4) 


where  62fc- 1  +  162*  =  8 -)£«*. 

The  least  squares  problem  is  to  find  b  such  that 

^(S£(x,y)-br*(x,y)]2 

D 

is  minimised,  and  where 

br  =  [a«, 6],  67, . . . ,  f>2;v] , 

*r  =  («'•..  ri.  V’li  -  •  - ,  V’2/v) , 


(3.5) 


b 


in  which 


02i-j(*.y)  =  e-J’,‘r  cos 27tj/  +  e~1',tX  cos2'>(ti/ 

^•2*(*,y)  =  ‘  (e-Jl*1’cos2')ty  -  e-2"'*1  cos¥jky) , 
and  /)  represent  the  rectangular  array  of  grid  points. 

The  least,  squares  problem  gives  rise  to  the  so  called  Siormal  equations" 


(3.6) 


Lb  =  r, 


(3.7) 


in  which 

Lij  =  ^2  C.(J-.yH>(.r,y), 

rr  =  ('V  r, .  i-> . r2lV].  (3.8) 

where 

r,  =  y^  S'  js,  y)c,(  j-.  u) 

The  discrete  least  squares  problem  suffers  from  the  loss  of  precision  due  to 
the  subtraction  of  two  similarly  valued  numbers  of  finite  word-length.  If  the  sam¬ 
pling  rate  is  sufficiently  high  or  grid  sparing  is  sufficiently  small,  then  all  sums 
may  be  replaced  by  integrals  with  an  error  proportional  to  the  square  of  the  grid 
spacing.  This  leads  to  a  procedure  known  as  Continuous  beast  Squares  where,  by 
performing  the  integrations,  it  may  be  shown  that  the  elements  of  the  symmetric 
matrix  L  are  given  by 

d-O’.')},  (39) 

1(1, 2fc  +  1)  =  ~/nt  |^p(l  -  e-2-”' )  J  , 

and  the  remaining  terms  are  more  easily  expressed  in  terms  of  the  supplemental 
function  F( 7m,7„)  such  that 

£(2m,  2n)  =  F(ym,  y„)  +  F(ym,%)  +  F(ym,  y„)  +  F( 7m,7„) 

and 

l(2m  +  l,2n+  1)  =  -F( 7m,7»)  +  f(7m,7„) +  F( 7m.7n)  -  ^(7m.7„). 

(3.10) 
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where,  for  w  4  », 


/’( Tin  ■■)„)  =  jj  >  2<lm+^”,x  ei»s'2iI„i/cos2')„i/<lr<ltf 

(i  _  ,-«■>- +-•••»')  r  i 

_  ii-  +-  i  -  +-  sin''''' +  '"! 

»\  /»*»  T  fn  }  l  fm  T  in 


“lui  "t"  "/  >»  )  "f"  '  */»»)]■ 


Note  that  for  in  =  /?, 


(  “nn  •  “lm  )  — 


(1  f si !)•>-, 


=  0  (since  sin  2-,n.  -f  2';j.  =  0). 


The  least  squares  matrix  for  the  problem  under  consider#) ion  becomes  ltifthly 
ill-comlit  ioned  as  tlio  mimUer  of  rows  (or  columns)  increases  beyond  7.  The  de¬ 
terioration  of  the  conditioning  of  L  is  illustrated  in  Table  I.  Such  behaviour 
means  that  the  inversion  of  L  on  a  computer  with  finite  word-lengths  may  not 
give  dependable  results  for  the  high  frequency  modes.  This  gives  some  idea  of  the 
maximum  number  of  parameters  which  can  be  reliably  estimated  in  the  inverse 
problem.  However,  because  these  parameters  can  be  strongly  correlated,  it  may 
be  necessary  to  include  higher  frequency  eigenfunctions  in  the  least  squares  model 
to  avoid  bias  in  estimating  the  low  frequency  components.  To  systematically  de¬ 
termine  the  order  of  the  model,  a  consistent  order  estimation  scheme  was  adopted 
(see  Ref.  [9])  which  involved  the  minimization  of  the  following  objective  function: 


Z  —  N  lno^  +  m  In  A, 


(3.13) 


where  N  is  the  total  number  of  data  mints,  m  is  the  number  of  parameters  assumed 
in  the  model,  and  ffm  ^  the  variance  of  the  least  squares  fit.  It  may  be  shown  that 
as  N  — »  oo,  m  may  be  determined  exactly. 

As  m  increases  beyond  7,  the  least  squares  matrix  may  deviate  from  being 
strictly  positive  definite  due  to  the  poor  conditioning.  To  overcome  this,  a  standard 
technique  which  involves  the  addition  of  a  small  number  c  to  the  diagonal  elements 
of  L  was  adopted  (see  Ref.  {10]).  Because  a  Choleski  decomposition  was  used  in 
the  matrix  solver,  the  value  of  t  should  be  of  the  order  of  the  square-root  of  the 
computer  truncation  errors.  A  value  of  t  =  10-6  was  selected  to  be  appropriate 
in  the  current  example. 


Number  of  Vs  12  3  1 

Number  uf  rows  in  L  3  5  7  9 

Condition  number^  II  SO  520  1070 

f  Thr  roinlition  iiihiiIkt  Inn1  is  I  l.li  m.  I  as  I  In’  of  llto  I.II  a's  I 

and  smallest  eigenvalue  of  (lie  tnalrix. 

Table  I.  Conditioning  of  the  least  s<|iiares  matrix. 


4.  RESULTS  AND  DISCUSSION 

The  minimization  of  the  bjeetive  function  expressed  by  e<|n  (3.13)  allows  a  sys¬ 
tematic  means  of  determining  the  number  of  parameters  which  should  be  included 
in  the  model.  The  procedure  adopted  was  to  compute  7  for  increasing  number  of 
Vs,  and  select  the  case  which  corresponds  to  a  minimum  7, .  For  the  two  noise  lev¬ 
els  considered  (±2.5  and  ±5  units),  this  procedure  returned  an  estimate  for  »n  of 
29  (or  14  Vs).  This  is  in  close  agreement  with  the  15  Vs  actually  used  to  generate 
the  simulated  SPATE  data.  However,  even  though  such  a  high  order  system  was 
considered,  most  of  the  parameters  were  estimated  with  substantial  errors.  Table 
2  shows  the  comparison  of  the  estimated  parameters  and  the  actual  values  used 
in  the  direct  problem  for  the  maximum-noisi  case.  It  is  seen  that  only  the  first  8 
parameters  were  estimated  reasonably,  which  is  consistent  with  the  conditioning 
of  the  least  squares  matrix  as  discussed  earlier  (see  Table  1).  Despite  this,  these 
parameters  were  able  to  be  reproduce  the  bulk  stress  field  with  exceptional  accu¬ 
racy  as  seen  in  Fig.  2  which  shows  a  comparison  of  the  bulk  stress  field  produced 
in  the  direct  problem  prior  to  the  addition  of  random  noise,  and  that  obtained  in 
the  inverse  problem. 

The  comparisons  of  the  <rzx ,  < 7vy  and  <tt,  fields  are  shown  in  Figs  3  to  6 
respectively.  Only  the  results  for  the  ±5  unit  noise  case  are  presented  herp  as  this 
scheme  was  found  to  be  relatively  insensitive  to  the  random  noise  present,  and  the 
the  results  for  the  ±2.5  unit  noise  case  were  essentially  the  same.  Interestingly,  it 
is  seen  that  the  inaccurately  predicted  higher  order  parameters  tend  to  affect  the 
quality  of  estimated  individual  stress  components  much  more  than  that  for  the 
bulk  stress  field.  The  reason  for  this  may  be  easily  seen  from  the  algebraic  form 
of  these  stress  fields.  For  example,  setting  x  —  0  in  eqn  (2.G),  and  comparing  this 
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with  o<(ii  (3.3),  if  can  l>o  seen  fliaf  (lie  stress  component  <rjr  has  terms  which  are 
proportional  to  *ij!  and  whilst  S  contains  only  terms  in  -lk  ami  Yk.  Therefore,  as 
k  (and  hence  y*  and  ^k )  hocomes  larger,  errors  in  <tk  would  become  more  and  more 
visible  on  stres-  components.  I’ort  iniately  though,  the  higher  order  parameters 
are  associated  with  terms  which  decay  rapidly  in  the  .r-dirccl ion.  and  as  shown 
in  these  figures,  the  spurious  stresses  are  routined  to  only  a  small  region  on  the 
loaded  edge,  lor  the  bulk  of  the  region  under  consideration,  the  inverse  scheme 
was  able  to  determine  all  three  components  of  stresses  accurately. 


k 

bk  ( Direct ) 

bk  (Inverse) 

1 

33.3 

33.3 

2 

48.7 

48.5 

3 

-4.13 

-3  96 

1 

35.2 

31.7 

r» 

-3  70 

-2  1 1 

6 

7.07 

7.00 

T 

-6.18 

-4.17 

8 

-11.9 

-12.1 

9 

-2.71 

-0.386 

10 

-12.8 

-9.62 

11 

1.88 

2.77 

12 

-3.12 

0.904 

13 

3.93 

-4.41 

14 

6.86 

1.37 

15 

2.31 

-5.59 

16 

8.07 

3.69 

17 

-1.26 

4.32 

18 

2.04 

5.16 

19 

-3.01 

0.733 

20 

-4.88 

-5.81 

21 

-1.96 

-0.742 

22 

-5.93 

0.167 

23 

1.00 

3.34 

24 

-1.55 

8.13 

25 

2.42 

-13.3 

26 

3.84 

-2,11 

27 

1.83 

-16.9 

28 

4.74 

-2.74 

29 

-1.13 

-3.62 

30 

0.748 

- 

31 

-1.30 

- 

Table  2.  Comparison  of  the  parameters  used  in  generating  S 
and  those  obtained  in  the  inverse  problem. 
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5.  CONCLUSION 

W  hilst  isopachic  data  such  as  those  obtained  in  a  thermoelast ic  stress  analysis 
alone  are  not  sufficient  for  individual  stress  components  to  be  deduced,  ihe  corre¬ 
sponding  stress  field  cannot  be  chosen  arbitrarily.  This  is  because  a  valid  stress 
field  must  satisfy  the  conditions  of  equilibrium  as  well  as  the  imposed  boundary 
conditions.  It  iias  been  demonstrated  that  by  making  use  of  these  conditions,  a 
least  squares  scheme  can  be  devised  (o  identify  the  stress  components  in  a  two 
dimensional  body.  For  the  case  considered,  the  scheme  was  found  to  be  relatively 
insensitive  to  random  noise,  although  due  to  ill-conditioning  of  the  least  squares 
matrix,  the  high  frequency  components  of  stresses  were  not  able  to  be  resolved. 
However,  these  high  frequency  components  of  stresses  dissipate  extremely  quickly 
in  accordance  with  St.  Yenant's  principle,  and  the  stress  components  in  t  he  greater 
part  of  the  region  under  considerat  ion  can  be  determined  successfully. 

For  the  purpose  of  demonstrating  such  a  combined  experimental-analytical 
stress  analysis  technique,  the  chosen  example  was  a  rather  simple  one.  In  practice, 
for  a  component  of  arbitrary  geometry,  the  problem  becomes  much  more  difficult 
and  resorting  to  the  analytical  solution  may  be  impractical.  In  such  a  case,  there 
is  scope  for  adopting  this  least  squares  approach  in  a  finite  elements  formulation. 
The  viability  of  this  combined  experimental-finite  elements  technique  is  currently 
under  investigation. 
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